Matching Conditions in Atomistic-Continuum Modeling of 

Materials 

Weinan and Zhongyi Huang^ 
^Department of Mathematics and PACM, Princeton University and 
School of Mathematics, Peking University 
^PACM, Princeton University and 
Department of Mathematical Sciences, Tsinghua University 

February 1, 2008 

Abstract 

A new class of matching condition between the atomistic and continuum regions is pre- 
sented for the multi-scale modeling of crystals. They ensure the accurate passage of large scale 
information between the atomistic and continuum regions and at the same time minimize the 
reflection of phonons at the interface. These matching conditions can be made adaptive if 
we choose appropriate weight functions. Applications to dislocation dynamics and friction 
between two-dimensional atomically flat crystal surfaces are described. 

Traditionally two apparently separate approaches have been used to model a continuous medium. 
The first is the continuum theory, in the form of partial differential equations describing the conser- 
vation laws and constitutive relations. This approach has been impressively successful in a number 
of areas such as solid and fluid mechanics. It is very efflcient, simple and often involves very few 
material parameters. But it becomes inaccurate for problems in which the detailed atomistic pro- 
cesses affect the macroscopic behavior of the medium, or when the scale of the medium is small 
enough that the continuum approximation becomes questionable. Such situations are often found 
in studies of properties and defects of micro- or nano- systems and devices. The second approach 
is atomistic, aiming at knowing the detailed behavior of each individual atom using molecular 
dynamics or quantum mechanics. This approach can in principle accurately model the underlying 
physical processes. But it is often times prohibitively expensive. 

Recently an alternative approach has been explored that couples the atomistic and continuum 
approaches. |[ ^, ^, |[ ^, 0]. The main idea is to use atomistic modeling at places where 
the displacement field varies on an atomic scale, and the continuum approach elsewhere. Two 
representative examples of such a coupled approach are the quasi-continuum method ^ and 
the coarse-grained molecular dynamics |^. Both methods have been successfully applied to 
quasi-static examples, but extension to dynamic problems has not been straightforward. The 
main difficulty lies in the proper matching between the atomistic and continuum regions. Since 
the details of lattice vibrations, the phonons, which are an intrinsic part of the atomistic model, 
cannot be represented at the continuum level, conditions must be met that the phonons are not 
reflected at the atomistic-continuum interface. Since the atomistic region is expected to be a small 
part of the computational domain, violation of this condition quickly leads to local heating of 
the atomistic region and destroys the simulation result. In addition, the matching between the 
atomistic-continuum interface has to be such that large scale information is accurately transmitted 
in both directions. 

The main purpose of this paper is to introduce a new class of matching conditions between 
atomistic and continuum regions. These matching conditions have the property that they allow 
accurate passage of large scale (scales that are represented by the continuum model) information 
between the atomistic and continuum regions and no reflection of phonon energy to the atomistic 
region. These conditions can also be used in pure molecular dynamics simulations as the border 



For the sake of clarity we will first explain the main issues and ideas on a simple problem: a 
one-dimensional chain of particles coupled by springs: 



ilj — Uj^i — 2uj + Uj-i (1) 
The spring constant is set to be 1. After discretization in time, we have 

' - ""+1 - + -"-1 (2) 

where u" is the displacement of the j-th particle at time t — nAt. 

Equation (^) is supposed to be solved for all integers j. Now let us assume that we will 
truncate the computational domain and only compute u" for j > 0. At j = 0, we will impose a 
new boundary condition to make sure that the phonons arriving from j > are not reflected back 
at j = 0. 

The phonon spectrum for (|^) is obtained by looking for solutions of the type u" = e''^"'^'^*+^^^ . 
This gives us the relation 

1 • ■ c 

— — sm = sni f {6) 
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At j = 0, we replace (||) by 

< 51 "fcJ^r'' "0,0 = (4) 

k.j>0 

We would like to determine the coefficients {akj}. For the simple problem at hand, it is possible 
to obtain analytical formulas of {a-kj} such that the imposition of (Q) together with the solution of 
(||) for J > reproduces exactly the solution of ^ if it was solved for all integer values of j, i.e. an 
exact reflectionless boundary condition can be found. The details of this is given in Q . This exact 
boundary condition should be the same as the one found numerically in . It represents the exact 
Green's function for (|^) which is nonlocal. However, this procedure appears to be impractical for 
realistic models, particularly when the atomistic region moves with time which is the case that 
interests us. Therefore we will not pursue this direction here. 

A practical solution is to restrict to a finite number of terms and look for the coefficients 
{flfcj } that minimize reflection. In order to do this, let us look for solutions of the type 

u] = giini^^t+ji) + i?(^)e'("'^^*-J«) (5) 

where lo is given by (|^). R{^) is the reflection coefficient at wavenumber ^. Inserting into (|^), 
we obtain 

The optimal coefficients {a^.j} are obtained by 

min / Wmnm^d^ (7) 







subject to the constraint 

R{0) = 0,i?'(0) = 0, etc. (8) 

Here W{^) is a weight function, which is chosen to be W{£,) = 1 in the examples below. 

Condition (^) guarantees that large scale information is transmitted accurately, whereas @) 
guarantees that the total amount of reflection is minimized. This procedure offers a lot of flexibil- 
ity. For example, instead of |i?((^)pc?^, we can minimize the total reflection over certain carefully 
selected interval. Another possibility is to choose the weight function to be the (empirically com- 
puted) energy spectrum. The coefficients {akj} may then change in time to reflect the change of 
the nature of the small scales. In practice, we found it preferable to use ^ with some 

small 6, instead of in order to minimize the influence of = tt for which we always 

have R{tt) = 1. 



-B- two layers 
-©- three-layers 




0.5 1 1.5 2 2.5 



Figure 1: Reflection coefficients for (^ij) and ([l^). 



Let us loofc at a few examples. If in we only keep the terms involving oi^o a-nd then 
imposing the condition i?(0) = gives 



(1 - A^X"^ + At7i^ 



n-l 



(9) 



If instead we keep terms involving ao,i,ai.o and ai^i, we can then impose both R{0) — and 
i?'(0) = 0. This gives us 

<^<-' + \^^K-'~u-) (10) 

Conditions of the type and ( p^ ) are intimately related to the absorbing boundary conditions 
proposed and analyzed in for the computation of waves. These conditions perform well for 

low wavenumbers but are less satisfactory at high wavenumbers. 

To improve the performance at high wavenumbers let us consider a case that include terms with 
k < 2,j < 3 and minimize ^ (with 6 = 0.1257r) subject to the condition i?(0) = 0, 

the optimal coefficients can be easily found numerically and are given by 



1.95264069 -7.4207 x lO^^ -1.4903 x lO^^ 
-0.95405524 7.4904 x lO^^ 1.5621 x 10"^ 
If instead we only include terms such that fc < 3, j < 2, then 

/ 2.9524 1.5150 x lO^^ 
-2.9065 -3.0741 x IQ-^ 



(11) 



(12) 



V 0.9541 1.5624 x IQ-^ / 



The resulting reflection coefficients R are displayed in Figure 1. 

We have applied this method to a number of problems. As the simplest model that encompasses 
most of the issues in a coupled atomistic/continuum simulation, we consider the Frenkel-Kontorova 
model 



— Xj+i - 2xj + Xj-i - U'{xj) + f 



(13) 



where J7 is a periodic function with period 1 , / is an external forcing. The continuum limit of this 
equation is simply the Klein-Gordan equation 
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Figure 2: Comparison of the displacement and velocity profiles computed using the full atomistic 
and the atomistic/continuum models, with / — 0.04. The left column shows the results in the 
whole computational domain. The right column shows the details near the dislocation. 



where K = U"{0). We consider the case when there is a dislocation and study its dynamics under 
a constant applied forcing. We take U{x) = {x — [xVf where [x] is the integer part of x. In this 
example we take (|l|) as our atomistic model, and ( |l4|) as our continuum model. For the coupled 
atomistic-continuum method, we use a standard second order finite difference method for ( [T^ ) 
in the region away from the dislocation, and we use (13) in the region around the dislocation. 
However, we also place finite difference grid points in the atomistic region. At these points, the 
values are obtained through averaging the values from the atomistic model. At the interface 
between the atomistic and continuum regions, we decompose the displacement into a large scale 
and a small scale part. The large scale part is computed on the finite difference grid, using (p^). 
The small scale part is computed using the reflectionless boundary conditions described earlier. 
The interfacial position between the MD and continuum regions is moved adaptively according to 
an analysis of the wavelet coefficients or the local stress. Both strategies lead to similar results. 
Care has to be exercised in order to restrict the size of the atomistic region. For example, when 
wavelet coefficients are used in the criteria to move the atomistic region, we found it more efficient 
to use the intermediate levels of the wavelet coefficients rather than the finest level. 

We first consider the case when a sharp transition is made between the atomistic and continuum 
regions with a 1:16 ratio for the size of the grids. Figure 2 is a comparison of the displacement 
and velocity fields computed using the full atomistic model and the coupled atomistic/continuum 
model, with / = 0.04. The atomistic region has 32 atoms. The full atomistic simulation has 10''. 
Dislocation appears as a kink in the displacement field. Notice that at the atomistic/continuum 
interface, there is still substantial phonon energy which is then suppressed by the refiectionless 
boundary condition. No reflection of phonons back to the atomistic region is observed. 

We next consider a case with / ~ 0.02, which alone is too weak to move the dislocation, but to 
the left of the dislocation, we add a sinusoidal wave to the initial data. The dislocation moves as a 
consequence of the combined effect of the force and the interaction with the wave. Yet in this case 
the same atomistic/continuum method predicts an incorrect position for the dislocation, as shown 
in Figure 3. The discrepancy seems to grow linearly in time. Improving the matching conditions 
does not seem to lead to significant improvement. 

The difference between this case and the case shown in Figure 2 is that there is substantially 
more energy at the intermediate scales. This is clearly shown in the energy spectrum that we 
comouted for the two cases but it can also be seen in Fieure 3 where an aDDreciable amount of 
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Figure 3: Comparison of the displacement and velocity profiles computed using the full atomistic 
and the atomistic/continuum models, with / = 0.02. The left column shows the results when the 
transition from the atomistic to continuum regions is sharp. The right columns shows the results 
when the transition is gradual. Solid line is the result of the atomistic/continuum method. The 
dash line is the result of the full atomistic method. Only the region near the dislocation is shown. 



a method that uses a sharp transition between the atomistic and continuum regions. We therefore 
consider the next alternative in which the atomistic/continuum transition is made gradually in a 
1:2 or 1:4 ratio between neighboring grids. The right column in Figure 3 shows the results of such 
a method that uses a gradual 1:2 transition. We see that the correct dislocation position is now 
recovered. 

Our second example is the friction between atomically flat crystal surfaces. To model this 
process atomistically, we use standard molecular dynamics with the Lennard-Jone potential 1 11 , 12 1. 
The two crystals are separated by a horizontal interface. The atoms in the bottom crystal are 
assumed to be much heavier than the atoms on top. A constant shear stress is applied at the top 
surface. The main issue here is how dissipation takes place. Physically the mean kinetic energy is 
dissipated through conversion into phonons which then convert into heat and exit the system. A 
standard practice in modeling such a process is to add a friction term to the molecular dynamics 



in order to control the temperature of the system |11, n2|. In contrast, we ensure the proper 



dissipation of phonons to the environment by imposing the reflectionless boundary conditions for 
the phonons. As a result we obtain a linear relationship between the mean displacement of the 
atoms in the top crystal as function of time, see Figure 4. The temperature of the system also 
saturates (Figure 4). Also plotted in Figure 4 is the result of the mean displacement computed 
using the combined atomistic/continuum method. Here the continuum model is the linear elastic 
wave equation with Lame coefficients computed from the Lennard- Jones potential. The agreement 
between the full atomistic and the atomistic/continuum simulation is quite satisfactory. 

In conclusion, we presented a new strategy for the matching condition at the atomistic/continuum 
interface in multiscale modeling of crystals. These conditions are adaptive if we choose the weight 
functions in (Q) to reflect the evolving nature of the small scales. They minimize the reflection of 
phonons and at the same time ensure accurate passage of large scale information. 

This work is supported by NSF through a PECASE award and by ONR grant N00014-01-1- 
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